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ABSTRACT 

Propagation  of  electromagnetic  energy  through  the  atmosphere  is  a  difficult  task 
because  of  temperature  fluctuations  and  index-of-refraction  inhomogeneities  which  de- 
grade the  beam's  coherence.  Understanding  this  phenomena  is  of  practical  importance 
for  optical  systems. 

This  thesis  presents  an  analytical  numerical  technique  which  simulates  the  effects 
of  atmospheric  turbulence.  The  extended  Huygens-Fresnel  principle  was  used  to  simu- 
late wave  propagation  in  a  two-dimensional  randomly  varying  medium,  which  is  repres- 
ented by  thin,  filtered.  Gaussian  phase  screens.  The  wave  optics  code  implements  both 
Fresnel  and  Fraunhofer  propagation,  by  employing  the  fast  Fourier  transform  (ITT) 
algorithm.  The  analytical  spatial  coherence  length,  p0  ,  and  normalized  intensity  vari- 
ance. G-jl\  of  the  perturbed  electric  field,  were  examined.  Results  support  the  concept 
of  intensity  saturation  for  weak  scattering  cases,  however,  differences  in  the  values  of  the 
theoretical  and  analytical  spatial  coherence  lengths,  occurred. 
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I.     INTRODUCTION 

Atmospheric  turbulence  degrades  the  coherence  of  electromagnetic  energy  propa- 
gating through  the  atmosphere.  The  refractive-index  inhomogeneities  associated  with 
the  turbulent  atmosphere  induce  phase  and  intensity  perturbations  across  the  wavefront. 
Understanding  the  propagation  and  scattering  of  optical  waves  in  random  media,  is  es- 
sential for  atmospheric  laser  beam  propagation  and  imaging  systems. 

This  thesis  models  the  propagation  through  a  random  media  by  means  of  the  ex- 
tended Huygens-Fresnel  principle.  A  two-dimensional,  thin,  Gaussian  phase  screen  re- 
presents the  randomly  varying  medium.  This  wave  optics  propagation  code  employs  the 
fast  Fourier  transform  (FFT)  algorithm  in  both  Fresnel  and  Fraunhofer  forms  of  prop- 
agation. The  Fresnel  region  incorporates  two  forms  of  the  diffraction  integral,  the 
transfer  function  form  for  "near  field"  distances,  and  the  convolution  form  for  "far  field" 
regions. 

Previously,  numerically  intensive  simulations  of  this  type  required  large,  super- 
computer systems.  However,  these  computations  were  performed  on  a  compact, 
desktop  computer  system. 

The  model's  accuracy  was  assessed  by  comparing  computer  values  of  the  spatial 
coherence  length.  p0  .  with  values  obtained  from  the  analytical  mutual  coherence  func- 
tion (MCF).  Additionally,  the  concept  of  intensity  variance  saturation  was  examined  for 
a  single  phase  screen  in  the  Fraunhofer  approximation.  The  normalized  intensity  vari- 
ance, a]jl2,  approaches  a  saturation  value  asymptotically  for  increasing  values  of  the 
index-of-refraction  structure  parameter.  On  .  Theoretical  calculations  suggest  saturation 
for  multiple  scattering,  however,  they  say  nothing  about  a  single  phase  screen  realiza- 
tion. This  thesis  provides  results  which  support  saturation  effects  for  single  scattering 
cases. 


II.     BACKGROUND 

A.     STATISTICAL  DESCRIPTION  OF  TURBULENCE 

1.  Random  Variables 

Maintaining  the  coherence  of  an  electromagnetic  wave  propagating  through  the 
atmosphere,  requires  an  understanding  of  the  effects  imposed  on  the  wave  by  the  tur- 
bulent medium.  The  fundamental  characteristic  of  atmospheric  turbulence  is  its  ran- 
domness, which  must  be  described  statistically.  In  addition  to  statistical  quantities, 
further  assumptions  must  be  made  of  atmospheric  turbulence.  These  include  the  con- 
cepts of  stationarity.  homogeneity  and  isotropy.  Stationarity  implies  that  random 
processes  are  time  independent.  Homogeneity  assumes  invariance  under  a  Galilean 
transformation  of  coordinates,  while  isotropic  variables  are  invariant  with  respect  to 
coordinate  rotations.  [Ref.  1]  Mathematically,  these  two  assumptions  imply  that  the 
statistics  at  two  points  r,  and  r2  depend  only  on  the  difference,  ru  =  \rx  —  r2\  • 

2.  Local  Homogeneity  and  Isotropy 

In  general,  atmospheric  random  variables  do  not  obey  the  assumptions  of 
stationarity.  homogeneity,  and  isotropy.  However.  Tatarski  [Ref.  2]  introduces  the 
concept  oi'  "local"  homogeneous  and  isotropic  random  variables.  This  concept  requires 
homogeneity  and  isotropy  within  a  localized  region  of  size  L0.  the  outer  scale  length. 
Furthermore,  the  difficulties  associated  with  nonstationary  random  variables,  are  re- 
moved by  considering  random  fields  with  stationary  first  increments  [Ref.  2].  Tatarski's 
method  applies  to  a  nonstationary  random  function  whose  mean  varies  slowly  with  time, 
by  considering  the  difference  of  the  function  at  two  different  locations.  The  slow  func- 
tional changes  do  not  affect  the  value  of  this  difference. 

3.  Structure  Functions 

Tatarski  introduces  the  structure  function 

DJ(?2-F])  =  <[!\72)-J{rl)f>,  (1) 

a  tensor  that  is  the  difference  between  two  quantities.  Some  important  aspects  of  the 
structure  function  arc;  that  its  general  form  is  valid  for  any  variable,  and  <  >  denotes 
an  ensemble  average  taken  over  all  possible  point  pairs  r2  .  rL  Assuming  homogeneity 
and  isotropy.  the  vector  dependence  reduces  to  the  magnitude  r12  =  \r2  —  r,|  ,  and  the 

structure  function  becomes 


DJ(r)  =  <[f[r2)-j{r])r>.  (2) 

Kolmogorov  showed  through  dimensional  analysis,  a  simple  power  law  dependence  of 
equation  (2).  over  a  limited  interval  called  the  inertial  subrange,  as 

Dj{r)  =  Cjr212.  (3) 

Tatarski  introduces  the  concept  of  "passive  additives",  which  are  quantities  independent 
of  position  in  the  vector  field,  and  do  not  directly  influence  the  dynamics  of  the  turbulent 
medium.   Temperature  is  a  passive  additive  and  has  a  structure  function  of  the  form 

Djir)  =  CJr2'3.  (A) 

As  long  as  r  remains  within  the  inertia!  subrange,  temperature  is  approximated  as  a 
passive  additive,  and  equation  (4)  is  valid.  Likewise  the  index-of-refraction  is  a  passive 
additive  with  a  structure  function 

Dn{r)=Cy!3.  (5) 

Since  the  index  of  refraction  depends  on  the  density  of  the  atmosphere,  Dr  and  DT  are 
related  by 

Dn  =  D7{l9x  lO-6-^- V.  (6) 

4.     Covariance,  Power  Spectral  Densities 

In  addition  to  the  structure  function,  other  characteristics  of  random  processes 
include  the  covariance  (or  correlation),  and  Power  Spectral  Densities.  It  is  the  interre- 
lationship of  these  three  quantities  which  provide  a  useful  method  for  analyzing  random 
processes.   The  covariance  between  two  random  variables  S  and  T  can  be  expressed  as 

Vst=  <T(r])  -  <T(ri)>][S(r2)  -  <S(r2)>]>.  (7) 

However,  more  frequently  it  is  the  autocovariance  function 

BTt=  <7Ir1)  -  <T(r,)>}\T(r2)  -  <T(r2)>]>,  (8) 

which  is  needed.     Furthermore,  if  T  is  homogeneous  and  <7>  =  0  is  assumed,  then 

equation  (S)  simplifies  to 


B1{r)  =  <T{rx)T{r1)>.  (9) 

Combining  this  relation  and  equation  (2),  gives  an  expression  for  the  relationship  be- 
tween the  structure  function  and  covariance  function  as 

D7{r)  =  2[Br7{0)-BTJ{y)).  (10) 

In  one  dimension  the  covariance  function  and  the  power  spectral  density  are  transform 

pairs  given  by 


IV(k) 


— (Kr 


e~lKrB{r)dr,  (11) 


and 


B(r)  =  j^  eiKr)Y{K)dK.  (12) 

"'-oo 

Using  the  fact  that  B(r)  is  an  even  function  and  substituting  equation  (12)  into  equation 
(10),  D\r)  becomes. 

D(r)  =  2  |       [1  -  cos{Kr)}U\K)clK.  (13) 

Tatarski  develops  an  expression  for  the  one-dimensional  Kolmogorov  spectral  density, 
given  by 

n-(K-)  =  0.1224Q;K"5/3.  (14) 

Discussion  to  this  point  has  been  of  one-dimensional  random  processes,  how- 
ever these  concepts  are  applicable  to  three-dimensional  cases.  Analogous  to  equation 
(11),  Tatarski  defines  the  three-dimensional  power  spectral  density  as 


/*+oo 


<D(ic) 


(*+oo     f+oo 


£  B(r)dr, 


— oo        — oo        — oo 


and  similarly,  the  correlation  function  is 


(15; 


f+oo     /*+oo    /*+oo 


B{7)  = 


(2«)J 


e^'F<l>(K)/K. 


-oo        — oo        —  oc 


Using  the  relation 


0)(k-: 


</ir 


2/TK  i/K 


the  three-dimensional  Kolmogorov  power  spectral  density  becomes 

(D(K-)  =  0.033  Qk~11/3. 


(16) 


(17) 


(18) 


B.     EM  PROPAGATION  THROUGH  TURBULENCE 
1.     Wave  Equation 

Based  upon  Tatarski,  Clifford  [Rcf.  3]  develops  theoretical  results  of  line- of- sight 
propagation  through  the  atmosphere,  directly  from  Maxwell's  equations.  Assuming  zero 
conductivity,  and  unit  magnetic  permeability  in  the  atmosphere,  as  well  as.  a  sinusoidal 
time  dependent  electromagnetic  field,  Maxwell's  four  equations,  in  Gaussian  units,  take 
the  form 


V  .  H  =  0, 


VxE  =  ikH, 


VxII  =  -iknE, 


(19.) 
(20) 
(21) 


V  .  (nzE)  =  0.  (22) 

Taking  the  curl  of  equation  (20),  and  substituting  it  into  equation  (21),  gives 

-V2£  +  V(\~.£)  =  A-V£.  (23) 


Rewriting  equation  (22)  in  the  form 

E  .  V/?2  +  /;2V  .  E  =  0,  (24) 

and  substituting  it  into  equation  (23).  yields  the  vector  form  of  the  wave  equation 

V2£  +  A-y£  +  2V(£.Vlog*)  =  0.  (25) 

The  third  term  of  equation  (25)  describes  the  change  in  polarization  of  a  propagating 
electromagnetic  wave.  This  term  is  negligible  as  long  as  the  wavelength  is  small  com- 
pared to  the  refractive  inhomogeneities.   Thus  equation  (25)  reduces  to 

V2E  +  k2n2E  =  0.  (26) 

Equation  (26)  is  the  vector  form  of  the  wave  equation  describing  propagation  through 
the  turbulent  atmosphere.  The  difficulty  in  solving  this  equation  lies  in  the  second  term 
containing  the  random  variable  n.  Various  methods  are  available  for  obtaining  solutions 
to  equation  (26),  each  of  which  relics  on  several  critical  approximations.  Strohbehn 
[Ref.  4]  lists  these  approximations  as: 

1.  Negligible  depolarization  efTects. 

2.  Negligible  back-scattering. 

3.  Use  of  the  parabolic  approximation  to  the  wave  equation. 

4.  Turbulence  is  uncorrelated  in  the  direction  of  propagation. 

2.     The  Method  of  Small  Perturbations-Born  Approximation 

Both  Tatarski  [Ref.  2]  and  Clifford  [Ref.  3].  solve  the  wave  equation  in  a  turbu- 
lent atmosphere  using  the  method  of  small  perturbations,  which  is  equivalent  to  the 
Born  approximation.  This  method  expands  the  electric  field  into  a  series  of  decreasing 
amplitudes,  and  the  refractive  index  into  a  power  series  in  the  form 

E=E0  +  E]  +  ....  (27) 

k  =  1  +  «j  + (28) 

Substituting  these  into  equation  (26)  and  equating  same  order  terms,  results  in  two 
equations 

V:£u  +  A2£0  =  0,  (29) 


V2£j  +  k2E]  +  2k2n]  £0  =  0,  (30) 

where  terms  of  order  n\  and  higher  are  ignored.  Assuming,  as  Tatarski  does,  that  the 
unperturbed  field  is  a  plane  wave  propagating  in  the  z-direction  represented  as 
£0  =  exp[ikz]  ,  equation  (30)  becomes 

V2E^  +  k2Ex  =  -2k2nxeik\  (31) 

an  inhomogeneous  partial  diilerential  equation  with  constant  coefficients.  Its  solution 
is  the  convolution  of  a  plane  wave  Green's  function  with  the  source  term,  or  inhomo- 
geneous term,  given  by 


£i,r>  =  17 


J 


,      exp(/A|;-  —  r'\)        0 
dr'         .-_^|         [2**«i(r#)  exp(ikz')}.  (32) 

s 


3.     The  Method  of  Smooth  Perturbations-Rytov  Approximation 

In  addition  to  the  Born  approximation,  Tatarski  develops  the  Rytov  approxi- 
mation which  assumes  a  solution  to  equation  (30)  of  the  form 

£=exp(xF)  =  exp(X  +  /S),  (33) 

or  simply 

E=Aexp{iS),  (34) 

where  A  is  the  amplitude  given  by  A  =  exp  X  .  Applying  equation  (33)  to  equation  (30) 
and  dividing  by  E0,  yields  the  Rytov  solution  given  by. 

V2£+  k2n\r)  =  V2  log  £+  (V  log  £)2  +  k2n2(r).  (35) 

Tatarski  further  shows  that  both  methods  of  approximation  are  equivalent. 

C.    HUYGENS-FRESNEL  THEORY 

As  we  have  seen,  the  theories  of  Tatarski  and  Clifford  use  the  diilerential  equation 
approach  to  solve  the  problem  of  propagation  through  a  turbulent  medium.  However, 
Lutomirski  and  Yura  [Ref.  5],  approach  this  problem  in  terms  of  integral  equations 
which  use  an  extended  Huygens-Fresnel  theory.  This  technique  is  equivalent  to  a  dif- 
ferential equation  approach,  but  it  is  easier  to  integrate  and  simulate  using  FFT  tech- 


niques.  Lutoinirski  and  Yura,  develop  an  extended  Muygens-Fresnel  theory  by 
introducing  a  random  phase  term  for  turbulence  in  the  Huygens-Fresnel  integral  which 
is  developed  in  standard  optics  texts  like  [Ref.  6].  This  additional  phase  perturbation 
takes  the  form  of  the  Rytov  approximation,  e*  .  The  extended  Huygens-Fresnel  integral 
is. 


£(F) 


-ik 


e(ik\r  —  r'\ 


\r  -  r 


t(r  )e        a  r 


(36) 


In  the  geometrical  optics  limit,  Fermat's  Principle  is.  Y-tkfn^zyiz.  With  a  power  series 
expansion  of  t,T .  equation  (36)  reduces  to  the  Green's  function  solution  of  Tatarski  and 
Clifford  given  in  equation  (32). 

Recently.  Martin  and  Flatte  [Ref.  7]  presented  an  atmospheric  turbulence  algorithm 
which  uses  the  differential  equation  approach  that  has  as  a  solution,  the  extended 
Huygens-Fresnel  integral.  A  filtered  random  Gaussian  phase  screen  introduces  the 
phase  perturbations  while  the  algorithm's  path  integral  method,  incorporates  a  multi- 
screen transfer  function  form  of  Fresnel  propagation. 

D.     MUTUAL  COHERENCE  FUNCTION 

The  effects  of  atmospheric  turbulence  can  be  expressed  in  terms  of  two  functions, 
the  Mutual  Coherence  Function  (MCF)  and  the  Modulation  Transfer  Function  (MTF). 
Lutomirski  and  Yura  [Ref.  5]  derive  the  first  concept  by  considering  the  average  intensity 
<!{}■)>  =  <E(T\)E[r1)>  .  of  equation  (36).  What  results  is  an  average  intensity  which  is 
a  product  of  the  autocovariancc  of  the  aperture  and  the  atmospheric  MCF.  The  atmo- 
spheric MCF  term  has  the  Rytov  form  <  exp(vI'"  -f  vr"')>  where  the  H1"  refers  to  the 
complex  phase  factor  at  the  r,  coordinate  and  VF  '  '  corresponds  to  the  r2  coordinate.  This 
term  is  log-normally  distributed,  as  long  as,  4'  is  composed  of  Gaussian  variables.  Using 
this  fact  and  results  in  Fried's  work  [Ref.  S  ].  the  atmospheric  MCF  was  written  in  terms 
of  the  wave  structure  function  D(p  ).  given  by 


<ey  '>  =  exp 
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where  D{p)  =  Dx(p)  +  D $([>).  Lutomirski  and  Yura  apply  the  structure  function  for  a 
plane  wave, 
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to  equation  (37)  which  can  be  written  as, 


<eCl"  +xv"]>  =  MCF{p)  =  exp 


P     \5/3 

Po 


(39) 


where  p0  ,  the  lateral  cohernece  length,  given  by 
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where  k  =  -^^  and  Q  is  the  index  of  refraction  structure  parameter  along  the  optical 
path  length  L.  represents  the  distance  where  the  spatial  coherence  of  the  wave  drops  to 
e~l  point  of  the  MCF.  [Ref.  8] 

E.     MODULATION  TRANSFER  FUNCTION 

Lutomirski  and  Yura's  concept  of  MCF  is  closely  related  to  the  MTF  of  Fried's 
[Ref.  S].  The  MCF  and  MTF  are  in  fact  the  same  function  but  expressed  in  terms  of 
different  variables.  The  MCF  is  measured  in  the  coordinates  of  the  propagation  field 
and  has  dimensions  of  distance,  while  the  MTF.  is  measured  in  the  image  plane  and  has 
the  dimensions  of  spatial  frequency.  Both  are  equivalent  under  a  transformation  be- 
tween the  two  planes  by  letting  p-*/.Rf  where  R  is  the  focal  length  of  the  optics  and  f  is 
the  spatial  frequency.  This  transformation  is  valid  under  the  Wiencr-Khintchine  theo- 
rem since  the  lens  Fourier  transforms  the  incident  electric  field  at  the  aperture  to  the 
image  plane.    Fried's  expression  for  the  atmospheric  long  term  MTF  is 


MTF{f)  =  exp 
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where  rB  =  2.1p0. 

In  calculating  the  MTF,  two  distinct  cases  exist.  These  are  a  short  term  and  long 
term  MTF.  The  short  term  exposure  describes  the  evaluation  of  the  wavefront  in  suffi- 
ciently short  time  intervals  such  that  the  turbulence  appears  fro/en.     The  long  term 


MTF,  is  a  single  long  time  integrated  exposure,  taking  into  account  every  turbulence 
configuration.  This  thesis  focuses  on  the  method  prescribed  by  the  short  term  MTF. 
[Ref  S] 

The  analytic  distinction  between  the  two  cases  lies  in  the  manner  in  which  the 
wavefront  distortions  are  handled.  Specifically,  the  distortion  attributed  to  a  random  tilt 
of  the  wavefront.  Tilt  results  from  the  varying  phase  fluctuations  across  the  aperture 
which  accumulate  along  the  optical  propagation  path.  These  fluctuations  produce  im- 
age motion  in  the  focal  plane  of  the  receiver.  For  a  very  short  exposure,  the  tilt  is  ex- 
tracted, by  fitting  a  mean  square  two-dimensional  flat  plane  to  the  electric  field  and 
rotating  it  through  an  angle  so  that  the  mean  wavefront  is  normal  to  the  direction  of 
propagation.  This  introduces  a  phase  shift,  resulting  in  the  displacement  of  each  curve 
about  the  optical  axis.  Fried's  [Ref.  S]  development  of  this  theory  suggests  a  higher 
MTF  at  all  spatial  frequencies  for  the  short  term  MTF.  The  short-exposure  MTF  for 
near-field  and  far-field  cases  is  respectively  given  by 
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where  t0  is  the  MTF  of  a  diffraction-limited  lens,  and  the  exponential  term  corresponds 
to  the  atmospheric  MTF.  These  equations  predict  a  near  field  short  term  MTF  that 
start  at  one.  declines  to  a  minimum,  then  increases  to  unity,  at  the  optical  cut  olf  fre- 
quency.   Figure  1  illustrates  this  phenomena. 

F.     DIFFRACTION  INTEGRAL 

This  section  presents  the  analytical  work  concerned  with  the  development  of  my 
propagation  algorithms.  It  begins  by  illustrating  the  approximations  that  were  used  to 
manipulate  the  diffraction  integral  and  then  proceeds  with  the  analysis  required  to 
transform  the  solution  of  the  Helmholtz  equation  into  two  forms  of  the  Huygens-Fresnel 
principle.  The  two  forms  are,  first,  a  convolution  form  suited  for  long  distance  propa- 
gation and  second,  as  a  transfer  function  form  for  short  distance  propagation.  Some 
initial  assumptions  made  by  Roberts  [Ref.  9J  include  the  following: 

1.  Light  propagates  in  the  k  direction. 

2.  The  wave  amplitude  is  known  in  the  xy-plane  at  z  =  0  . 
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Figure  1.        Short  Term  Mutual  Coherence  Function  for  an  8  x  8  Subaperture. 
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3.  Polarization  effects  are  negligible. 

4.  The  electric  field  amplitude  is  a  scalar  function  V(r,z). 

Following  Roberts,  the  analysis  begins  with  the  Huygens-Fresnel  integral  for  the 
propagation  of  light  waves  given  as 


r(F'z)  =  77 


</pK(p,0)exp 


r  2m       '  2  ,   ,-      -,2l 
xp|^-j-V2   +I'-Pl   J' 


(43) 


where  A  is  the  wavelength  of  light  and  p  is  a  vector  in  the  aperture  plane,  r  is  a  vector 
in  the  image  plane,  and  z  is  the  propagation  direction.  From  the  Fresnel  approximation 
where,  r  <  <  z  and  p  <  <  z,  factor  z2  from  the  square  root  and  expand  by  a  binomial 
expansion.  Equation  (43)  becomes 
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Since  the  exponential  phase  factor  outside  of  the  integral  does  not  affect  intensity 
measurements,  equation  (44)  is 


l'(r.z)^—. —    dpV{p, 0)  exp 


in   1-      —  ,2 

77  ir-pl 


(45) 


Expanding  the  quadratic  term  gives, 


V(r,  z)  =  -yi-  f^rip,  0)  exP[  -f-  [r2  -  2F .  p  +  p2] 

u    J  \_    s.Z 

=  — - -  exp    -7— r       dpVCp,  0)  exp    — —  p      exp 
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This   is   the   convolution   form   of  the   Fresnel   integral,  which  is   equivalent   to   the 
Fraunhofer  integral  except  for  the  quadratic  phase  factor  in  the  integral. 

To  obtain  the  transfer  function  form,  the  analysis  begins  with  denoting  V{f,z)  as  the 
Fourier  transform  of  V{r,z)  given  by 
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V(f\z)  =        drexp(-27tif.r)V{r,z), 


(47) 


where  V(r,z)  is  given  in  equation  (45).    Interchanging  the  order  of  integration  yields 


V(f,z)  = 


—i 
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Next,  a  chanse  of  variables  is  made  where 


r   =  r  —  p. 


Substituting  this  into  equation  (4S)  gives  the  following  equation 


(49) 


V(f,z)  = 
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The  r  '  integral  is  replaced  with  its  Gaussian  transform  pair, 

/'/.z  exp(  -inAzf), 
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The  inverse  Fourier  transform  of  equation  (52)  yields 


V{7,  z)  =  jdfexp{2nif.r)  exP(  -in).zf)^dpV(p,  0)  exp(  -2mf-p).  (53) 

This  expression  is  the  transfer  function  form  of  the  diffraction  integral,  which  is  equiv- 
alent to  the  solution  of  the  differential  equation  approach  used  by  Martin  and  Flatte 
[Ref.  7J.  Reviewing  the  form  of  equations  (46)  and  (53).  the  need  for  two  equations  de- 
scribing different  propagations,  is  obvious.  In  one  instance,  the  propagation  distance  z 
enters  in  the  denominator  of  the  exponential  term.  It  is  at  long  distances  that  the  ex- 
ponential term  varies  slowly.     On  the  other  hand,  equation  (53)  is  suited  for  short 
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propagation  distances,  as  z  enters  into  the  numerator  of  the  exponential  for  slow  vari- 
ations in  this  case. 
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III.     NUMERICAL  SIMULATION  MODEL 

This  numerical  simulation,  modeled  wave  optics  propagation  of  an  electromagnetic 
wave  through  a  random  medium,  represented  by  a  two-  dimensional  Gaussian  phase 
screens.  Techniques  in  this  model  required  certain  "tools"  and  their  testing.  These 
"tools"  included  a  need  for  a  reliable  random  number  generator,  used  in  generating  the 
random  phase  screens,  and  an  efficient  two-dimensional  FFT,  used  in  approximating  the 
diffraction  integrals.    But  first,  a  discussion  of  the  experimental  arrangement  is  needed. 

A.  EXPERIMENTAL  ARRANGEMENT 

Due  to  the  extensive  numerical  calculations  in  this  simulation,  a  Compaq  deskpro 
80386-20  computer  was  used.  It  features  a  64  megabyte  hard  drive  and  16  megabytes 
of  memory.  In  addition  to  the  20-MHz  80387  coprocessor,  a  Weitek  1167  math 
coprocessor  was  added  to  enhance  execution  speed.  The  32  bit  Fortran-386  compiler 
was  from  Silicon  Valley  Software  (SYS)  and  uses  Phar  Lap  Software  to  extend  the  op- 
erating system  beyond  one  megabyte.  The  math  and  graphics  packages  were  produced 
by  Scitech  Scientific.  The  Compaq  has  a  640  x  480  pixel  YGA  graphics  monitor.  Both 
the  HP  Laser  Jet  Scries  II  and  Panasonic  KX-P1092i  multi-mode  printers  were  used  in 
this  arrangement. 

B.  COMPUTER  PRELIMINARIES 

1.     Random  number  generator 

The  main  purpose  of  a  computer  simulation  is  to  approximate  natural  phe- 
nomena. To  make  things  realistic,  random  number  sequences  were  used  to  introduce 
stochastic  variations.  One  might  ask,  what  minimal  criteria  should  a  particular  random 
number  generator  satisfy?  Certainly  the  most  important  criterion  is  that  the  sequence 
of  numbers  is  sufficiently  random.  Other  criteria  are  uniformity,  reproducibility,  mini- 
mum memory,  fast,  non-repeating,  and  statistically  independent. 

The  more  difficult  characteristic  to  satisfy  is  statistical  independence.  Thus  a 
series  of  tests  were  needed  to  provide  a  quantitative  measure  of  the  generator's  per- 
formance. There  are  two  kinds  of  statistical  tests:  empirical  tests  and  theoretical  tests. 
Empirical  tests  focus  on  how  the  computer  manipulates  groups  of  numbers  from  the 
sequence  and  evaluates  certain  statistical  quantities.    Perhaps  the  best  known  of  all  sta- 
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tistical  tests  are  the  "Chi-Square"  tests.    Theoretical  tests,  on  the  other  hand,  establish 
characteristics  of  a  sequence  using  methods  based  on  recurrence  rules.  [Ref.  10] 

For  the  purpose  of  this  simulation,  only  empirical  tests  were  applied  to  the 
SVS-Fortran  random  number  generator  called  Ran^Ij  .  The  following  five  tests  were 
considered: 

1.  Frequency  Test.    This  test  determines  whether  or  not  the  sequence  of  numbers  are 
uniformly  distributed  as  U(0,1). 

2.  Serial  Test.    This  test  is  an  extension  of  the  frequency  test  to  two  dimensions  or 
matrix  form. 

3.  Lagged-Product  Test.  This  test  checks  for  correlations  between  successive  numbers 
over  a  given  lag  period. 

4.  Run  Tests: 

a.  Runs  up  and  down.    This  tests  for  long  increasing  and  decreasing  sequences  of 
numbers. 

b.  Runs  above  and  below  the  mean.    This  tests  for  long  sequences  with  values 
consistently  above  or  below  the  mean. 

The  results  of  the  five  tests  arc  provided  in  Appendix  A.  Of  these  five  tests,  the 
Lagged-Product  and  Runs  Tests  are  the  most  critical  when  simulating  atmospheric  tur- 
bulence. These  three  tests  determine  whether  or  not  correlation  is  introduced  from  the 
random  number  generator  which  produces  erroneous  results  in  the  simulation.  The 
SVS-Fortran  random  number  generator  met  the  test  criteria  and  proved  to  be  one  of  the 
better  generators.  However,  this  generator  has  one  significant  draw  back.  The  random 
sequences  begin  at  one  of  two  different  values  depending  on  whether  the  seed  value  is 
positive  or  negatiave.  Thus  it  is  critical  that  the  random  numbers  be  called  continuously 
in  a  loop  to  avoid  restarting  the  sequence,  thereby  introducing  unwanted  correlation. 
2.     Fast  Fourier  Transform 

The  most  repeatedly  used  algorithm  throughout  the  numerical  simulation  was 
the  fast  fourier  transform.  Therefore,  it  was  necessary  to  use  the  most  efficient  algorithm 
available.  The  Scitech  Scientific  math  package  provided  several  options,  with  subroutine 
FFT2C  .  best  suited  for  this  numerical  simulation.  This  subroutine  uses  a  complex  array 
input.  The  other  FFT  considered  was  a  routine  coded  by  Dr.  Walters  which  he  received 
in  a  demonstration  package  provided  by  Infotek.  This  FFT  utilizes  real  arrays  and  will 
henceforth  be  refered  to  as  subroutine  FFT  . 

Each  subroutine  was  timed  for  various  configurations.  Specific  timing  results 
are  contained  in  Appendix  B.     Some  general  results,  however,  are  that  subroutine 
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FFT2C  was  faster  for  one-dimensional  cases,  with  subroutine  FFT  faster  for  the  two- 
dimensional  cases.  The  decline  in  efficiency  of  FFT2C  can  be  attributed  to  the  extra 
coding  required  to  convert  between  real  arrays  and  complex  arrays.  Subroutine  FFT 
was  selected  over  FFF2C  since  the  simulation  utilized  the  two-dimensional  form  of  the 
FFT. 

Other  techniques  which  were  employed  to  increase  the  efficiency  included  in- 
stalling a  Weitek  coprocessor  in  the  Compaq.  This  reduced  the  processing  time  to  ap- 
proximately one-third  that  of  the  original  time.  The  use  of  common  blocks  vice 
dimension  statements  further  decreased  processing  time  by  5%.  Finally,  a  portion  of  the 
FFT  algorithm  was  modified  from  wim  =  s'm(ang)  ,  to  wim  =  x  1  —  wre2  ,  with  a  negli- 
gible decline  in  efficiency  by  0.01%. 

Because  of  the  discrete  nature  of  the  simulation,  there  exists  problems  and  lim- 
itations associated  with  implementating  the  FFT.  One  such  problem  was  that  of  clas- 
sical edge  diffraction  associated  with  the  phase  or  amplitude  discontinuities  at  the  edges 
of  the  finite  screen.  As  the  propagation  distance  z  increases,  the  edge  diffraction  spreads 
toward  the  center  of  the  screen  making  more  and  more  of  the  diffraction  pattern  erro- 
neous. Buckley  [Ref.  1 1]  defines  the  distance  from  the  ends  of  the  screens  where  the  edge 
diffraction  is  important  as 

D(z)  =  2V^7  +  200*2,  (54) 

where  z  is  the  propagation  distance  and  </>„  is  the  root  mean  square  phase  deviation  im- 
posed on  the  wave  by  the  screen.  The  severity  of  edge  effects,  however,  is  reduced  by 
the  aliasing  introduced  in  the  FFT  implementation.  Aliasing  transforms  the  linear 
screen  to  a  "circular"  one  with  the  last  point  associated  with  the  first.  This  results  in  a 
continuity  of  phase  and  amplitude  at  the  edges  of  the  screen  [Ref.  11]. 

The  most  important  limitation  was  the  finite  spatial  range  imposed  by  the 
maximum  available  grid  size.  This  places  a  constraint  on  the  available  range  of  fre- 
quencies used  in  the  FFT  from  the  lowest  given  by  f!uin=\jL  ,  to  the  highest, 
./max  —  "/-£.  where  L  is  the  grid  length  and  n  is  the  number  of  grid  points.  Figure  3  on 
page  20  illustrates  this  setup.  The  FFT  provides  a  least  squares  fit  of  sine  and  cosine 
fuctions  to  the  perturbed  wavefront  phases.  However,  this  method  prevents  an  accurate 
representation  of  the  wavefront  at  low  frequencies  for  a  Kolmogorov  k~11/3  power  spec- 
trum. To  alleviate  this  problem,  a  subaperture  was  superimposed  at  the  center  of  the 
grid.    The  subaperture  helps  low  frequencies,  which  contain  a  large  portion  of  the  am- 
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plitude,  but  hurts  the  high  frequencies  by  limiting  the  inner  scale.  There  exists  some  high 
frequency  edge  effects,  however,  these  are  minimal. 

C.     PROCEDURE 

Figure  2  provides  a  summary  of  the  coding  contained  in  Appendix  C.  This  con- 
ceptual diagram  illustrates  an  overview  of  the  procedural  steps  of  the  Fraunhofer  prop- 
agation algorithm. 

1.  Input  Parameters 

Since  this  step  is  straightforward,  extensive  discussion  is  not  required.  However, 
it  is  important  to  note  that  the  input  parameters  were  both  fixed  and  variable.  The  array 
size  and  Filter  value  were  fixed  quantities,  but  the  subaperture  size,  seed  value,  On  value, 
and  propagation  distance,  took  on  different  values.  The  actual  variable  names  are  doc- 
umented at  the  beginning  of  the  code. 

2.  Aperture  Mutual  Coherence  Function 

The  second  step  in  this  procedure  calculates  the  aperture  MCF.  Subroutine 
MCF  does  this.  In  this  subroutine  the  initial  wavefront  was  represented  in  the  computer 
as  an  L  x  L  square  array  of  complex  numbers.  Centered  within  this  complex  electric 
field  was  a  subaperture.  The  initial  complex  electric  field  had  a  value  of  zero,  every- 
where, except  for  the  real  part  of  the  subaperture,  which  had  the  value  of  one. 
Figure  3  illustrates  this.  With  the  electric  field  created,  it  was  direct  Fourier  transformed 
(DFT)  by  subroutine  DFTIFT  .  The  intensity  of  the  electric  field  was  calculated,  and 
then  inverse  Fourier  transformed  (I FT),  yielding  the  aperture  MCF. 

3.  Planar  Electric  Field 

The  complex  electric  field  was  created  by  the  same  method  prescribed  in  sub- 
routine MCF  .  It  is  important  to  realize  that  the  concept  of  aperture  si/e  was  used  in 
two  ways.  One  way  corresponds  to  the  simulated  aperture  while  the  other  pertains  to 
an  aperture  with  physical  dimensions.  The  simulated  aperture  is  actually  a  matrix  or 
grid  which  directly  reflects  the  dimensioning  size.  For  example,  the  simulated  L  x  L 
complex  electric  field  was  actually  a  two-dimensional  matrix  corresponding  to  a  256  x 
256  two-dimensional  array.  From  the  input  parameters,  the  simulated  subaperture  takes 
on  grid  sizes  ranging  from  an  8  x  8  to  200  x  200  square  matrix. 

The  second  referencing  to  an  aperture  refers  to  an  aperture  with  actual  physical 
units.  The  physical  subaperture  was  assigned  a  value  of  0.3125  meters,  which  remains 
fixed  regardless  of  the  simulated  subaperture  size.    The  length  L  of  the  complex  electric 
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Figure  3.        Three  Dimensional  Perspective  of  the  Initial  Electric  Field. 
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field,  on  the  other  hand,  had  variable  physical  lengths  depending  on  the  simulated 
supaperture  size.   The  physical  value  of  L  was  calculated  from 

L  =  m  x  — : —  [meters],  (55) 

isize 

where  m  is  the  physical  subaperture  length  in  meters,  nr  is  the  integer  value  of  the  array 
dimensioning,  and  isize  is  the  integer  value  for  the  simulated  subaperture. 
4.     Phase  Screens 

a.  Generation 

Phase  screens  were  created  in  subroutine  GGAUS  by  constructing  a  L  x  L 
matrix  of  Gaussian  distributed  random  numbers.  Each  position  was  assigned  an  inde- 
pendent random  number  n.  thereby,  requiring  n2  random  numbers.  Since  the 
SVS- Fortran  random  number  generator  was  uniformly  distributed,  an  algorithm  pro- 
vided by  Knuth  [Ref.  10]  transformed  the  distribution  into  a  Gaussian  one.  Two  inde- 
pendent, two-dimensional  real  arrays  called  phaser  and  phase!  were  created,  which 
represent  the  real  and  imaginary  components  of  a  two-dimensional  complex  phase 
screen.    In  this  algorithm,  the  imaginary  part  of  the  phase  screen  was  set  to  zero. 

The  domain  in  which  the  phase  screens  are  created  is  arbitrary,  however, 
filtering  was  done  in  the  Fourier  plane.  Creating  the  phase  screen  in  frequency  space 
vice  real  space,  reduces  the  requirement  for  an  additional  FFT  when  transforming  from 
real  space  to  frequency  space.  Martin  and  Flatte  [Ref.  7]  proposed  this  method,  but  it 
creates  difficulties  in  absolute  normalization.  Further  discussion  is  presented  in  Chapter 
four.  This  simulation,  on  the  other  hand,  generated  the  random  Gaussian  phase  screen 
in  real  space.  This  in  turn  was  DFT'd  to  frequency  space  where  the  complex  phase 
screen  was  filtered  and  then  the  filtered  phase  screen  was  IFT'd. 

b.  Filtering 

Phase  screens  were  filtered  to  obtain  the  correct  power  law  form.  The  fil- 
tering function  used  in  subroutine  FLTR  was, 

O0(k)  =  2nk25x<\\,  (56) 

where  6X  is  the  slab  thichness  and  <1>„  =  0.033QK-"3  ,  which  gives  the  relationship  be- 
tween the  phase  spectrum  and  refractive-index  spectrum  [Ref.  2.  pp.  101-102.].  Filtering 
was  accomplished  by  multiplying  each  phase  screen  spectrum  by  the  square  root  of 
equation  (56).  The  correct  filtering  method  requires  circular  frequency  filtering  instead 
of   a    linear    one,    because    of   two-dimensional    isotropy.       The    correct    form    is 
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k  =  N  K;  +  k\  ,  which  is  radial  everywhere  except  at  the  origin,  where  it  is  zero  [Ref.  12]. 
This  was  reflected  in  the  simulation  by  setting  the  position  (1,1)  equal  to  zero  in  each 
two-dimensional  real  array  which  makes  up  the  complex  phase  screen. 

It  is  important  to  realize  that  although  the  filtering  concept  is  simple  and 
straightforward,  the  actual  implementation  is  not.  The  difficulty  arises  from  the  sym- 
metry properties  of  the  digitally  filtered  phase  screen. 

A  one-dimensional  case  offers  a  simple  illustration  of  this  concept.  The 
FFT  of  a  complex  function  which  contains  only  real  components,  results  in  a  symmetric 
function  in  frequency  space  about  the  Nyquist  frequency  for  the  real  components,  and 
an  anti-symmetric  function  for  the  imaginaries.  With  these  symmetries  present  in  the 
frequency  domain,  the  correct  implementation  of  the  filtering  is  to  mimick  these  sym- 
metries. Thus  a  "folding"  technique  about  the  Nyquist  frequency  was  required.  How- 
ever, when  this  concept  was  extended  to  a  two-dimensional  case,  as  in  the  phase  screen 
in  the  simulation,  the  symmetries  imposed  by  the  FFT.  were  no  longer  apparent  in  the 
frequency  domain.  To  simplify  the  symmetry  requirements,  a  real  phase  screen  was 
used.  The  real  and  imaginary  spectral  components  were  filtered  by  folding  about  the 
Nyquist  frequency.  It  was  suggested  by  both  Brigham  [Ref.  13]  and  Martin  and  Flatte 
that  a  complex  phase  screen  containing  both  real  and  imaginary  components,  yields  two 
entirely  distinct  phase  screens.  However,  nothing  was  provided  to  support  this  hypoth- 
esis. 

Another  important  consequence  of  filtering  resides  in  the  units.  The  phase 
screens  were  filtered  in  k  units,  but  the  ITT  algorithm  operates  in  frequency  units. 
Therefore,  it  was  necessary  to  make  a  change  of  variables  prior  to  applying  the  I  IT. 
The  relationship  used  for  the  change  of  variables  is  k  =  27:v  . 

The  paper  by  Martin  and  Flatte  [Ref.  7],  specifies  an  additional  normaliza- 
tion factor  of  A"1  in  equation  (56),  where  Ak  =  -rr-  and  where  N  is  the  number  of  grid 
points  and  A  is  the  sampling  interval.    Martin  and  Flatte  provided  no  explanation  for 
the  additional  A;1  in  the  filtering.  It  was  not  included  in  the  filtering  code. 
5.     Implementation 

After  the  filtered  phase  screen  was  IFT'd  into  real  space,  it  was  introduced  into 
the  code  as  a  phase  screen  and  multiplied  with  the  electric  field.  The  array  phaser.  which 
contains  the  desired  phase  field,  takes  the  form  of  the  Rytov  approximation,  e16  .  in  the 
extended  Huygens-Fresenel  integral.  This  form  assumes  that  only  phase  perturbations 
and  not  amplitude  variations,  occur. 
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6.     Propagation  Methods 

As  indicated  in  previous  sections,  the  Huygens-Fresnel  technique  was  used  to 
simulate  the  propagation  of  light.  This  was  accomplished  by  applying  the  FFT  to  the 
perturbed  electric  field.  Both  the  "far-field"  and  "near-field"  propagation  methods  were 
considered. 

a.  Fraunhofer 

Of  the  two  different  propagation  methods,  the  single  screen  Fraunhofer 
propagation  is  by  far  the  simplest  technique  to  implement  and  the  one  implemented  in 
this  thesis.  The  uniform  coherent  plane  wave  at  the  aperture  was  FFT'd  yielding  the 
desired  diffraction  pattern.  Looking  more  closely,  one  can  see  that  under  certain  cir- 
cumstances. Fraunhofer  propagation  is  just  a  special  case  of  the  long  distance  convo- 
lution form  of  Fresnel  propagation  given  by  equation  (46).  There  exists  two  situations 
when  this  occurs.  One  is  a  "far-field"  case  for  large  distances,  where  the  point  of  obser- 
vation is  at  infinity.  The  other  case,  is  when  a  spherical  curvature  is  placed  on  the  wave 
at  the  aperture.  This  curvature  cancels  the  quadratic  phase  factor  in  the  Fresnel  form, 
at  the  focal  point. 

Both  Fresnel  and  Fraunhofer  algorithms  were  needed  for  propagation.  It 
was  essential  to  verify  and  validate  each  case  before  building  on  the  pre-existing  codes. 
The  Fraunhofer  algorithm  provided  the  basis  of  this  simulation.  Since  the  Fraunhofer 
diifraction  pattern  is  well-known,  it  provided  a  means  to  verify  the  existing  code  by 
comparing  the  simulated  diffraction  pattern  with  theoretical  results.  Figure  4  and  Fig- 
ure 5  illustrate  one  three-dimensional  quadrant  of  the  diffraction  pattern  o^  an  unper- 
turbed electric  field,  for  two  different  subaperture  sizes.  While  Figure  6  illustrates  one 
three-dimensional  quadrant  of  a  perturbed  electric  field  difiraction  pattern  for  a  16  x  16 
subaperture. 

b.  Fresnel 

Although  the  Fresnel  propagation  forms  were  implemented  but  not  tested 
in  this  thesis,  discussion  is  warranted  since  multi-screen  Fresnel  propagation  is  predom- 
inantly used  in  thermal  blooming  and  multiple  scattering  scenarios.  Propagation 
through  the  turbulent  boundary  layer  also  requires  Fresnel  propagation  codes.  Two 
different  forms  are  used  for  Fresnel  propagation,  which  apply  a  straightforward  FFT  to 
evaluate  them.  These  two  forms  however,  do  not  allow  for  a  variable  receiving  array 
size.  In  addition,  choosing  the  correct  number  of  sample  points  is  essential.  An  effective 
approach  used  to  resolve  these  problems  is  the  Fresnel  number. 
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Figure  4.         Three  Dimensional  Diffraction  Pattern  for  an  8  .\  8  Subaperture. 
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Figure  5.        Three  Dimensional  Diffraction  Pattern  of  a  16  x  16  Subaperture. 


25 


<*A 


Figure  6.         Perturbed  Diffraction  Pattern  of  a  16  x  16  Subaperture. 
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The  Fresnel  number  is 

Hf-^f*-.  (57) 

where  A  r  is  the  receiving  aperture  size,  A  p  is  the  transmitting  aperture  size.  /  is  the 
wavelength  of  light,  and  z  is  the  propagation  distance.  Implementation  of  the  correct 
Fresnel  form  depends  on  the  Fresnel  number.  When  the  Fresnel  number  is  smaller  than 
the  number  of  grid  points  in  the  field  length,  that  is,  Nf<N  ,  the  long  distance  propa- 
gation algorithm  is  used.  Conversely,  when  N>Nf  ,  the  short  distance  code  is 
implemented. [Ref.    14  ] 

The  long  distance  propagation  code  uses  the  convolution  form  of  the 
diffraction  integral.  Implementation  requires  placing  a  curvature  on  the  electric  field 
wavefront  at  the  aperture.  The  subroutine  called  quad!  ,  does  this.  Multiplying  the 
phase  screen  with  the  electric  field,  gives  a  perturbed  electric  field  that  propagates  the 
entire  distance  by  means  of  one  FFT.  In  the  Fourier  plane,  the  quadratic  phase  factor 
called  quad2  scales  the  electric  field. 

The  transfer  function  form,  on  the  other  hand,  is  suited  for  short  distances. 
In  this  case,  the  entire  propagation  distance  is  divided  into  equally  spaced  slabs.  Two 
FFT's  are  required  to  propagate  the  distance  of  each  slab.  This  is  accomplished  by  the 
following  method: 

1.  Mesh  the  electric  field  and  phase  screen. 

2.  Apply  the  direct  FFT. 

3.  Multiply  the  field  by  the  propagation  transfer  function  subroutine  called  trnsfr  . 
A.  Apply  the  inverse  FFT. 

This  procedure  is  repeated  for  each  phase  screen  until  the  observing  plane  is  reached. 
[Ref  7] 

Both  forms  of  Fresnel  propagation  were  included  in  the  simulation  code 
contained  in  Appendix  C,  however,  neither  form  of  Fresnel  propagation  was  exercised. 
7.     Atmospheric  Mutual  Coherence  Function 

The  final  step  in  this  simulation  calculated  the  atmospheric  MCF.  The  same 
procedure  used  to  calculate  the  aperture  MCF,  was  applied  to  the  perturbed  electric 
field,  with  one  exception.  The  difference  is  that  division  of  the  composite  atmospheric- 
optical  MCF  by  the  aperture  MCF,  gives  the  atmospheric  MCF. 
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IV.     RESULTS 

The  main  objectives  of  this  thesis  were  to  demonstrate  that  the  simulation  gives  re- 
sults that  provide  insight  to  weak  scattering  wave  propagation  and  examine  the  accuracy 
of  the  simulation  code.  The  limitations  imposed  by  the  computer  mechanics  and  the 
actual  method  of  implementing  turbulence  theory  are  discussed.  Specific  areas  of  inter- 
est include  the  filter  function,  MCF,  and  saturation  effects. 

A.     FILTERING 

As  previously  mentioned  in  chapter  three,  Martin  and  Flatte  proposed  a  very  dif- 
ferent approach  to  creating  the  random  Gaussian  phase  screen.  Their  method  suggests 
creating  the  phase  screen  in  frequency  space,  vice  real  space,  to  reduce  the  required 
number  of  FFT's  from  two  to  one.  Since  the  domain  in  which  the  phase  screen  is  gen- 
erated, is  arbitrary,  this  approach  seems  plausible.  However,  this  method  proved  to  be 
awkward.  As  the  subaperture  size  was  successively  doubled,  it  was  necessary  to  increase 
the  strength  of  turbulence.  Q,  by  a  factor  of  ten  each  time,  in  order  to  produce  the 
identical  level  of  turbulence  as  in  the  previous  phase  screen.  Additionally,  since  the 
phase  screen  was  created  in  frequency  space,  and  no  symmetries  were  present,  the  fil- 
tering technique  did  not  reflect  any  folding  about  the  Nyquist  frequency.  It  is  not  clear 
that  aliasing  was  accounted  for  in  the  Martin  and  Flatte  algorithm.  Hence,  the  IFF  of 
the  phase  screen  appears  to  result  in  a  statistically  incorrect  phase  screen.  Additionally, 
it  is  not  clear  how.  with  one  FFT.  Martin  and  Flatte  handled  the  2r  and  1  N  normal- 
ization requirements.  With  a  round  trip  of  FFT's,  the  normalization  problems  are  au- 
tomatically handled.  The  weakly  filtered  phase  screen,  in  turn,  led  to  MCF  curves  which 
grossly  overestimated  the  coherence  length,  p0.  These  problems  obtained  from  Martin 
and  Flatte  '  s  approach  to  simulating  turbulence  led  to  the  current  coding  which  created 
one  phase  screen  in  real  space  and  applied  a  folding  about  the  Nyquist  frequency,  in  the 
filtering,  to  account  for  the  symmetries  introduced  by  the  FFT's. 

Kolmogorov  theory  of  atmospheric  turbulence  predicts  that  the  graph  of  the  phase 
screen  power  spectral  density  versus  k  yield  a  slope  of —11/3.  Figure  7  shows  that 
equation  (56)  produced  a  filtered  phase  screen  that  reflects  the  -11/3  slope,  as  well  as, 
the  correct  folding  technique.  Identical  slopes  were  expected  for  all  subapertures.  as  well 
as,  all  possible  angles  which  reflect  the  circular  filtering.  Other  subaperture  profiles  show 
a  consistent  slope  value  of -11/3.    Figure  8  corresponds  to  a  32  x  32  subaperture  at  a 
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45  degree  angle  while  Figure  9  is  representative  of  a  64  x  64  subaperture  at  a  90  degree 
angle.  Isotropy  is  apparent  in  that  the  circular  filtering  was  implemented  correctly 
within  the  error  introduced  by  the  randomness  in  the  screen  and  the  ability  to  linearly 
fit  a  line  through  the  data  points. 

B.     MCF 

The  MCF  of  the  electric  field  provides  one  method  of  analyzing  the  accuracy  of  the 
simulation.  To  verify  that  the  MCF  was  correctly  computed,  the  aperture  MCF  of  an 
unperturbed  electric  field  was  calculated.  This  was  easily  accomplished  since  the  image 
intensity  and  aperture  MCF  are  transform  pairs  and  are  analytical  for  simple  square  and 
circular  apertures.  The  aperture  MCF  is  just  the  autocorrelation  of  the  aperture  func- 
tion which  is  evaluated  by  calculating  the  area  of  overlap  of  two  identical  apertures  as 
they  are  moved  laterally  apart.  For  a  square  subaperture.  the  autocorrelation  yields  the 
triangle  function,  with  maximum  value  of  1.0  and  minimum  value  of  0.0  corresponding 
directly  to  the  subaperture  size.  Figure  10  illustrates  the  MCF.  or  autocorrelation,  of 
an  8  x  S  and  16  x  16  square  subaperture. 

The  MCF  corresponding  to  the  atmospheric  turbulence  was  determined  by  dividing 
the  MCF  of  the  perturbed  electric  field  by  the  aperture  MCF.  The  value  of  p0.  the  spa- 
tial coherence  length,  was  determined  from  the  turbulence  MCF  curve.  Figure  11  il- 
lustrates a  simulated  pu  value  of  3.20  mm  for  a  64  x  64  subaperture  with  On  —  LvlO-13. 
The  theoretical  value  calculated  from  equation  (40)  yields  p0  =  2.41mm.  Although  the 
64  x  64  subaperture  gave  accurate  results  other  subaperture  configurations  did  not. 
When  Q  increased,  the  MCF  curves  fell  oil  rapidly  towards  zero  for  all  subapertures. 
Figure  12  illustrates  this  for  a  64  x  64  subaperture.  All  subaperture  configurations  were 
run  for  two  dillerenct  Q  values.  The  simulated  p0  values  were  plotted  against  the  cor- 
responding subapertures  for  each  case.  Figure  13  reflects  the  Q  —  l.rl0~13p0  values,  and 
Figure  14  is  for  Q  —  IjcIO-14.  The  desired  trend  is  for  the  simulated  p0  values  to  approach 
the  theoretical  value  as  the  subaperture  size  increased.  This  trend  is  visible  in 
Figure  13  where  Q  =  LvlO-13  and  p0  =  2.4mm.  However.  Figure  14  shows  continuously 
decreasing  p0  values  past  the  theoretical  value  of  9.6mm.  The  results  of  figures  thirteen 
and  fourteen  point  to  a  problem  that  may  involve  edge  effects  or  aliasing,  resulting  from 
undersampling.  The  p0  values,  which  were  on  the  order  of  millimeters,  were  smaller  than 
the  tens  of  centimeter  distances  of  the  subaperture  mesh  size. 

In  an  attempt  to  pinpoint  the  problem,  the  simulation  code  was  changed  to  increase 
the  number  of  frequencies  in  the  subaperture  by  using  a  512  x  512  array.    The  results 
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Figure  7.        Filtered  Phase  Screen  of  a  16  x  16  Subaperture. 
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Figure  8.         Filtered  Phase  Screen  of  a  32  x  32  Subaperture. 
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Figure  9.         Filtered  Phase  Screen  of  a  64  x  64  Subaperture. 
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Computed  Autocorrelation  of  a  Square  Subaperture. 
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Figure  11.         Mutual  Coherence  Function  of  a  6-4  x  64  Subaperture. 
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Figure  12. 


Mutual  Coherence  Function  Trends. 
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Subaperture  Coherence  Lengths  for  Q  =  l.rlO-1J. 
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Figure  14.        Subaperture  Coherence  Lengths  for  Q  =  IatIO-'4. 
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were  identical  to  those  for  the  256  x  256  dimensioning,  within  the  arithmetic  error  of  the 
algorithm.  The  information  provided  by  the  512  dimensioning  was  that  the  inaccuracy 
of  the  code  was  not  due  to  an  inner  scale  problem,  however,  a  low  frequency,  outer  scale 
problem  may  still  exist. 

C.     SATURATION 

The  last  area  of  investigation  was  whether  or  not  the  simulation  predicts  the  satu- 
ration of  intensity  for  an  extended  medium  modeled  by  a  single  phase  screen  realization. 
Saturation  is  generally  considered  to  be  caused  by  multiple  scattering,  a  scattered  wave 
interfering  with  a  distorted  wave.  One  would  expect  from  the  Rytov  approximation  of 
turbulence  theory,  that  saturation  will  occur  even  in  Fraunhofer  propagation.  This  as- 
sumption stems  from  the  representation  of  turbulence  in  the  form.  e'e,  which  has  a 
magnitude  bounded  between  plus  and  minus  one.  Figure  15  illustrates  that  the  nor- 
malized intensity  variance  saturates  with  increasing  turbulence.  Theoretically,  a  nor- 
malized variance  of  one  is  expected  for  Rayleigh  statistics.  However,  it  is  premature  to 
assume  that  saturation  is  inherent  in  Fraunhofer  cases,  until  the  MCF  results  are  veri- 
fied. 
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Figure  15.         Intensity  Variance  Saturation  of  a  16  x  16  Subaperture. 
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V.     CONCLUSIONS  AND  RECOMMENDATIONS 

This  thesis  simulated  the  propagation  of  plane  waves  through  an  extended  two- 
dimensional  random  media  using  a  single  phase  screen  technique.  The  atmospheric 
turbulence  had  a  Kolmogorov  k-11'3  pure  power  law,  while  the  propagation  was  strictly 
Fraunhofer.  Limitations  of  its  applicability  were,  primarily,  the  finite  spatial  range  im- 
posed by  the  available  grid  size.  Diffraction  patterns,  correlation  functions  and  intensity 
variance  saturation  at  the  observation  plane,  were  investigated. 

The  results  provided  by  the  simulation  suggest  general  agreement  with  the  turbu- 
lence theory.  Saturation  for  weak  scattering  was  supported  by  this  model.  The  MCF 
curves,  although  not  completely  correct,  provided  insight  to  the  theory  and  illustrated 
problems  which  are  still  present  in  the  current  coding.  Some  specific  problems  include, 
potential  errors  in  the  implementation  of  the  filtering  technique  from  edge  effects,  alias- 
ing and  inner  scale  problems,  or  from  incorrect  normalization. 

Any  further  research  on  this  topic  should  begin  with  resolving  the  inaccuracies  still 
present  in  the  simulation  coding.  Several  possible  reasons  were  presented,  however,  an 
error  in  the  filtering  of  the  phase  screen  seems  to  be  the  most  likely  cause.  Further 
testing  can  be  conducted  on  the  phase  screen,  to  include  calculating  the  phase  screen 
variance,  as  well  as.  its  structure  function  De.  By  comparing  the  simulated  phase  screen 
with  the  theoretical  structure  function  for  /)„.  this  technique  will  verify  whether  or  not 
the  simulated  phase  screen  accurately  represents  turbulence. 

An  assumption  was  made,  that  aliasing  was  not  a  problem,  since  the  problems  as- 
sociated with  undcrsampling  were  not  apparent.  Aliasing  can  be  tested  by  using  finer 
grid  sizes  and  observing  changes  induced  by  the  higher  spatial  frequencies. 

After  the  coding  is  working  correctly,  both  the  convolution  and  transfer  function 
form  of  Fresnel  propagation  can  be  implemented  and  exercised.  In  addition,  the  array 
sizes  should  be  modified  to  incorporate  a  dimension  of  512  or  1024.  Finally,  during  the 
phase  screen  generation,  phasci  should  be  filled  with  random  numbers  to  give  two  usable 
phase  screens.  It  is  not  obvious  whether  two  independent  phase  screens  will  be 
produced,  or  whether  the  total  energy  will  be  distributed  among  the  two  phase  screens. 
Therefore,  testing  should  be  conducted  to  ensure  that  each  usable  phase  screen  possesses 
the  correct  statistics. 


40 


APPENDIX  A.     RANDOM  NUMBER  GENERATOR  TEST  DATA 

The  following  table  provides  the  statistical  results  of  the  five  empirical  tests  run  on 
the  random  number  generator,  RanrIj.  The  X2  values  were  determined  from  the  Chi- 
Square  table  in  Bevington  [Ref  15  ]. 

The  results  of  the  Lagged  Product  test  correspond  respectively  to  the  theoretical 
mean,  /ir  ,  calculated  mean,  \x  ,  theoretical  standard  deviation,  aT  ,  and  calculated 
standard  deviation,  o.   A  lag  of  three  was  tested. 

Table   1.      RANDOM  NUMBER  GENERATOR  TEST  DATA 


TEST 

DEGREES  OE 
EREEDOM 

X- 

PROBABILITY 

Frequency 

9 

4.4 

88.3°  o 

Serial 

20 

17.1 

65.0% 

Up  and  Down 

S 

4.1 

84.5" o 

Above  and  Below  the 
Mean 

12 

6.2 

90.5°  o 

TEST 

uT 

M 

a; 

a 

Lagged  Product 

0.250 

0.259 

0.100 

0.092 
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APPENDIX  B.     FFT  TIME  SERIES  DATA 

This  appendix  contains  the  results  of  the  comparison  between  the  two  FFT  sub- 
routines, FFT2C  and  subroutine  FFT  .  A  function  routine  called  SCNDS,  was  imple- 
mented in  the  code  to  provide  accurate  time  measurements. 

Table  2.     FFT  TIME  SERIES  DATA 


COMPUTER  SETUP 

FTT 

FFT2C 

1  DIMENSION 

TIME(sec) 

TIME(sec) 

2,;  with  20-MHz  80387 

2.64 

2.64 

21'  with  20-MHz  80387 

25.54 

25.26 

2r  with  20-MHz  80387 

54.05 

53.44 

2  DIMENSION 

TIME(sec) 

TIME(sec) 

128  x  12S  with  20-MHz  80387 

6.86 

7.08 

256  x  256  with  20-MHz  803S7 

31.30 

?2.0s 

256  x  256  with  Weitek 

12.08 

12.14 

256  x  256  with  Weitek  and  common 
block 

11.48 

11.56 

512  x  512  with  Weitek  and  common 

block 

49.54 

52.40 
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APPENDIX  C.     SIMULATION  CODE 

The  simulation  code  contained  in  this  appendix  incorporates  Fraunhofcr  and  both 
forms  of  Fresnel  propagation.    This  thesis  only  exercised  the  Fraunhofer  propagation. 


/-i  -'-  .'-  J.  Jm  J-  -J-  .J-  Jm  J*  Jm  - '-  ju  -•  -  - '  -  JL  JU  J-  J*  -r  -  J-  -'-  -'-  y  -  y-  -'-  -*  -  - '-  -'-  J-  -  ■-  -'-  - '-  *'-  -  *  -  -.*-  J  -  -'-  - '.  ..'-  **-  »U  -'-  -'-  -  •-  -' -  -•-  - '-  »*-  *'-  -v  y-  -* .  J*  -•-  -'-  -  ■-  -  •-  -'-  -'-  » -  -•  -  *i  -  y-  -'-  y  -  -f  - 

C 

C  THIS  CODE  PROVIDES  A  QUALITATIVE  VIEW  OF  BOTH  FRAUNHOFER  AND 

C  FRESNEL  DIFFRACTION  BY  OBSERVING  THE  PERTURBATION  IMPOSED  ON  A 

C  MONOCHROMATIC  PLANE  WAVE  PROPAGATING  THROUGH  A  TURDULENT 

C  MEDIUM.   THE  TURBULENT  MEDIUM  IS  INTRODUCED  IN  THE  FORM  OF  A 

C  STOCHASTIC  PHASE  SCREEN.   PROPAGATION  OF  THE  ELECTRIC  FIELD  IS 

C  ACCOMPLISHED  THROUGH  FFT'S. 

C 

C  *********MMr^^ 

c 

C  GLOSSARY  OF  VARIABLE  NAMES: 

C 

C  1.   ARRAYS: 

C 

C  RE  -  ONE  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR ,  WHICH  IS  USED  TO 

C  MANIPULATE  THE  REAL  PART  OF  THE  COMPLEX  ELECTRIC  FIELD  IN 

C  THE  FFT  SUBROUTINE.   THIS  ARRAY  IS  REPEATEDLY  USED  THROUGHOUT 

C  THE  CODE. 

C 

C  RIM-  ONE  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR ,  WHICH  IS  USED 

C  TO  MANIPULATE  THE  IMAGINARY  PART  OF  THE  COMPLEX  ELECTRIC 

C  FIELD  IN  THE  FFT  SUBROUTINE.   THIS  ARRAY  IS  REPEATEDLY  USED 

C  THROUGHOUT  THE  CODE. 

C 

C  FIELDR  -  TWO  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR 

C  CONTAINING  THE  REAL  PART  OF  THE  COMPLEX  ELECTRIC 

C  FIELD.   THIS  ARRAY  IS  REPEATEDLY  USED  THROUGHOUT  THE  CODE 

C 

C  FIELDI  -  TWO  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR 

C  CONTAINING  THE  IMAGINARY  PART  OF  THE  COMPLEX  ELECTRIC 

C  FIELD.   THIS  ARRAY  IS  REPEATEDLY  USED  THROUGHOUT  THE  CODE 

C 

C  FILL  -  TWO  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR  USED  AS  A 

C  DUMMY  ARRAY  IN  THE  I FT  OF  THE  POWER  SPECTRUM  WHICH  YIELDS 

C  THE  MCF. 

C 

C  FIELDM  -  TWO  DIMENSION  ARRAY  OF  DIMENSION  NR  X  NR  REPRESENTING 

C  THE  MAGNITUDE  OF  THE  PERTURBED  ELECTRIC  FIELD. 

C 

C  FMCF  -  TWO  DIMENSION  ARRAY  OF  DIMENSION  NR  X  NR  REPRESENTING  THE 

C  INTENSITY  OF  THE  PERTURBED  ELECTRIC  FIELD.   THIS  ARRAY  IS 

C  USED  IN  DETERMINING  THE  MCF. 

C 

C  FNORM  -  TWO  DIMENSION  ARRAY  OF  DIMENSION  NR  X  NR  REPRESENTING  THE 
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C  INTENSITY  OF  THE  UNPERTURBED  ELECTRIC  FIELD.   THIS  ARRAY 

C  IS  ALSO  USED  IN  DETERMINING  THE  MCF. 

C 

C     PHASER  -  TWO  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR  CONTAINING 

C  THE  REAL  PART  OF  THE  RANDOM  COMPLEX  PHASE  SCREEN. 

C 

C     PHASEI  -  TWO  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR  CONTAINING 

C  THE  IMAGINARY  PART  OF  THE  RANDOM  COMPLEX  PHASE  SCREEN. 

C 

C     FMAG  -  ONE  DIMENSION  SLICE  OF  THE  PERTURBED  MCF  USED  IN  THE 

C  GRAPHICS  ROUTINE.   REAL  ARRAY. 

C 

C     DIST  -  ONE  DIMENSION  REAL  ARRAY  REPRESENTING  THE  PIXELS 

C  CORRESPONDING  TO  THE  VALUES  IN  THE  ARRAY  FMAG. 

C 

C     VARAIBLES: 

C 

C     NR  -  DIMENSION  OF  THE  ARRAYS  EXACTLY  AS  SPECIFIED  IN  THE  DIMENSION 

C  STATEMENTS  IN  THE  CALLING  PROGRAM.   INPUT  INTEGER.  INDICE. 

C 

C     N2  -  ONE  HALF  OF  NR.  INTEGER. 

C 

C     M  -  POWER  OF  2  IN  THE  DIMENSIONING.   USED  IN  THE  FFT  SUBROUTINE. 

C 

C     ISIZE  -  INTEGER  VALUE  CORRESPONDING  TO  A  PARTICULAR  CHOICE  FOR  AN 

C  APERTURE  SIZE.  INPUT  VARIABLE. 

C 

C     NSIZE  -  INTEGER  VALUE  CORRESPONDING  TO  THE  ACTUAL  APERTURE  SIZE. 

C 

C     DELMSH  -  REAL  VALUE  REPRESENTING  THE  SAMPLING  INTERVAL  FOR  A 

C  PARTICULAR  APERTURE  SIZE. 

C 

C     SEED  -  REAL  INPUT  VARIABLE  USED  TO  BEGIN  A  RANDOM  SEQUENCE  OF 

C  NUMBERS  FOR  THE  SUBROUTINE  GGAUS. 

C 

C     NYES  -  INPUT  INTEGER  VARIABLE  WHICH  SELECTS  WHETHER  TURBULENCE  IS 

C  INTRODUCED  IN  THE  CODE. 

C 

C     FILTER  -  FIXED  INPUT  VALUE  USED  IN  THE  FILTERING  FUNCTION.   REAL. 

C 

C     CN2  -  REAL  INPUT  VARIABLE  REPRESENTING  THE  INDEX-OF-REFRACTION 

C  STRUCTURE  PARAMETER.   DETERMINES  THE  AMOUNT  OF  TURBULENCE 

C  INTRODUCED  IN  THE  FILTERING  FUNCTION. 

C 

C     DREC  -  REAL  INPUT  VARIABLE  REPRESENTING  THE  RECEIVING  FIELD  SIZE. 

C 

C     DTRNS  -  FIXED  VALUE  FOR  THE  TRANSMITTING  FIELD  SIZE.  REAL. 

C 

C     Z  -  REAL  INPUT  VARIABLE  REPRESENTING  THE  TOTAL  PROPAGATION 

C         DISTANCE  OF  THE  ELECTRIC  FIELD. 

C 

C     DELX  -  REAL  VARIABLE  FOR  THE  PROPAGATION  DISTANCE  TO  EACH  SLAB. 

C 

C     NUMSCR  -  INTEGER  INPUT  VARIABLE  USED  IN  NEAR-FIELD  FRESNEL 

C  PROPAGATION.   THIS  VARIABLE  CORRESPONDS  THE  THE  NUMBER 

C  OF  EQUALLY  SPACED  SLABS  WHICH  MAKE  UP  THE  TOTAL  DISTANCE. 
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c 

C     WVL  -  FIXED  VALUE  FOR  THE  WAVELENGTH. 

C 

C     ICNT  -  INTEGER  VALUE  WHICH  DETERMINES  WHETHER  OR  NOT  THE  FIELD  HAS 

C  PROPAGATED  THE  TOTAL  DISTANCE  Z.   USED  IN  NEAR-FIELD 

C  FRESNEL  PROPAGATION. 

C 

C     PI  -  VALUE  OF  PI. 

C 

C     TPI  -  TWICE  THE  VALUE  OF  PI. 

C 

C     MODE  -  REAL  VALUE  WHICH  DETERMINES  THE  FORM  OF  PROPAGATION. 

C 

C     SIGN  -  REAL  VALUE  EITHER  1.0  OR  -1.0  WHICH  DETERMINES  WHETHER  THE 

C  FFT  IS  DIRECT  OR  INDIRECT.   IT  ALSO  DETERMINES  WHETHER 

C  NORMALIZATION  OCCURS. 

C 

C     FLDM  -  MAXIMUM  VALUE  IN  THE  PERTURBED  MCF  ARRAY. 

C 

C     FMAX  -  MAXIMUM  VALUE  IN  THE  PERTURBED  ELECTRIC  FIELD  ARRAY. 

C 

C     GRAPHICS: 

C 

C     THE  FOLLOWING  VARIABLE  NAMES  ARE  EITHER  SPECIFIC  TO  THE  SVS 

C     GRAPHICS  ROUTINE  OR  USED  TO  MANIPULATE  DATA  FOR  GRAPHING: 

C 

C     NDEX,  IREG,  ANS,  GETC,  IUNITP,  IUNITV,  ISYMB,  ITNO,  MON,  NPRIN, 

C     MODE,  ONE,  TITLE,  IX,  IY,  XMAX ,  ICOLOR,  INCR,  NTOT. 

C 

C     *************************^^ 

C 
C 

COMMON  /BLK1/  RE( 256) ,RIM( 256) 

COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDI( 256 ,256) 

COMMON  /BLK3/  PHASER( 256 , 256) ,PHASEI( 256 , 256) 

DIMENSION  NDEX(20),FIELDM(256,256),FMCF(256,256),FNORM(256,256) 

DIMENSION  FMAG(130),DIST(130),FILL(256,256) 

INTEGER-2  IREG(9) 

DOUBLE  PRECISION  PI 

DATA  NDEX/ 15 ,15,7,7,8,8,14,14,5,5,4,4,2,2,3,3,1,1,0,0/ 

DATA  RE/256-0. 0/ ,RIM/256*0. 0/ 

DATA  FMAG/ 130*0. 0/ ,DIST/ 130*0. 0/ 

CHARACTERS  ONE 

CHARACTERS  1  TITLE 

CHARACTERS  ANS , GETC 

DATA  TITLE/' THE  SEED  VALUE  IS=  '/ 

DATA  ONE/'     '/ 

DATA  IUNITP/ 10/, IUNITV/20/ 

0PEN(4,FILE='LN.DAT' , STATUS= ' NEW ' ) 

ISYMB=22 

ITNO=5 

MON=18 

PI=3. 141592653589792 

NPRIN=0 

MODE=0 
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999  CONTINUE 

THIS  SECTION 
SIMULATION. 


') 


OF  THE  PROGRAM  SETS  UP  THE  INPUT  PARAMETERS  FOR  THE 


*)' HELLO. . 
*)' SEVERAL 

*)' 

*)'THE  VARIABLE 


*) 'FIRST  VALUE  TO  ENTER.   INPUT  THE  INTEGER  VALUE. 
)NR 

*)' 

*)'THE  SECOND  VARIABLE  OF  INTEREST  IS  "NSIZE" 
*)' VARIABLE  DIMENSIONS  THE  SIZE  OF  THE  PLANAR 
*)' FIELD.  SELECT  ONE  OF  THE  FOLLOWING.  ' 


*)* 
*)' 


1. 
2. 
3. 
4. 
5. 
6. 
7. 


WRITE (* 
WRITE (* 

WRITE(* 

WRITE (* 

WRITE (* 

WRITE (* 

READ(*, 

WRITE (* 

WRITE (* 

WRITE (* 

WRITE (* 

WRITE (* 

WRITE (* 

WRITE(* 

WRITE (* 

WRITE  (* 

WRITE (* 

WRITE (* 

WRITE (* 

READ(*,*)ISIZE 

IFCISIZE.EQ.  1)  THEN 

NSIZE=50 

DELMSH=.  0031 

INCR=3 
ELSEIF(ISIZE.EQ 

NSIZE=32 

DELMSH=.  0049 

INCR=1 
ELSEIF(ISIZE.EQ 

NSIZE=i6 

DELMSH=.  0098 

INCR=2 
ELSE1F(ISIZE.EQ.  4) 

NSIZE=8 

DELMSH=.  019531 

INCR=2 
ELSEIF(ISIZE.EQ. 5) 

NSIZE=4 

DELMSH=. 039063 

INCR=1 
ELSEIF(ISIZE.EQ. 6) 

NSIZE=2 

DELMSH=. 0781 

INCR=1 
ELSE 

NSIZE=1 

DELMSH=. 

INCR=1 
ENDIF 

VRITEOV 
WRITE  (*,' 


LET  US  BEGIN  THIS  SIMULATION  BY  ENTERING 
INPUT  PARAMETER  VALUES.  ' 


WHICH  DIMENSIONS  THE  ARRAY  SIZE  IS  THE 


THIS' 
ELECTRIC' 


FOR 
FOR 
FOR 
FOR 
FOR 
FOR 
FOR 


100  X  100 
64  X  64' 
32  X  32' 
16  X  16' 
8X8' 
4X4' 
2  X  2' 


2)  THEN 


3)  THEN 


THEN 


THEN 


THEN 
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ANOTHER  INPUT  PARAMETER  IS  THE  SEED  VALUE  OF  THE 
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WR I TE( *,*) 'RANDOM  NUMBER  GENERATOR.   INPUT  THE  SEED  VALUE  OF  +1.0* 
WRITEC*,*) 'OR  -1.  0' 
RE ADC*,*) SEED 

WRITEC*,*)' 

WRITE(*,*) 'FINALLY  INPUT  AN  INTEGER  OF  VALUE  1  FOR  TURBULENCE,  OR' 

WRITEC*, *)'0  FOR  NO  TURBULENCE.  ' 

READ(*,*)NYES 

IF(NYES.EQ.  1)  THEN 

WRITE(*,*)' 

WRITEC*,*) 'INPUT  THE  VALUE  OF  FILTER' 

READ(*,*)FILTER 

ENDIF 

WRITEC*,*)' 

WRITEC*,*) 'INPUT  THE  VALUE  FOR  CN2' 

READC*,*)CN2 

WRITEC*,*)'  ' 

c     WRITEC*,*) 'INPUT  THE  RECEIVING  FIELD  SIZE  IN  METERS' 
c     READC*,*)DREC 

WRITEC'--,*)' 

WRITEC*,*) 'INPUT  THE  TOTAL  PROPAGATION  DISTANCE  IN  METERS' 

READ(*,*)Z 

WRITEC*,*)'  ' 

c     WRITEC*,*) 'INPUT  THE  NUMBER  OF  SCREENS  TO  BE  USED  EITHER  FOR' 
c     WRITE(*,*) 'FRESNEL  OR  FRAUNHOFER  ' 
c     READC*,*)NUMSCR 

IF(SEED.  GE.  1.0)  THEN 
0NEC1: 2)='+l' 
ELSE 
0NEC1:  2)='-l* 

ENDIF 

TITLE(20: 21)=0NE(1: 2) 

N2=NR/2 

WVL=.  5E-6 

M=AL0GCREAL(NR))/AL0G(2. ) 

DTRNS=. 3125 

ICNT=0 
C 

C     THE  FOLLOWING  STATEMENT  DETERMINES  WHICH  FORM  OF  FRESNEL  IS  TO  BE 
C     USED 
C 

c     MODE=INT( ( 2*DTRNS*DREC ) / C  WVL*Z ) ) 
C 

C     THE  FOLLOWING  SUBROUTINE  CALLED  MCF  DETERMINES  THE  AUTOCORRELATION 
C     OF  THE  APERTURE  WHICH  IS  TO  BE  USED  IN  DETERMINING  THE  ATMOSPHERIC 
C     COHERENCE  LENGTH. 
C 

CALL  MCF(FNORM,NR,M,PI,NSIZE,N2,DELMSH) 
C 

C     THIS  SECTION  CREATES  THE  PLANAR  ELECTRIC  FIELD 
C 

DO  40  1=1, NR 
DO  40  J=1,NR 

FIELDRCI,J)=0.  0 
FIELDI(I,J)=0.0 
40  CONTINUE 

DO  41  I=N2-NSIZE+1,N2+NSIZE 
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DO  41  J=N2-NSIZE+1,N2+NSIZE 
FIELDR(I,J)=1.  0 
41  CONTINUE 
C 

C     THIS  SECTION  DETERMINES  THE  SLAB  THICKNESS(ES)  FOR  WHICH  THE 
C     ELECTRIC  FIELD  IS  PROPAGATED  THROUGH. THIS  IS  VALID  FOR  BOTH 
C     FORMS  OF  PROPAGATION. 
C 

C     IF(NR.  LE.MODE)  THEN 
C      DELX=Z/NUMSCR 
C      ELSE 

DELX=Z 
C 

C     THE  FOLLOWING  SUBROUTINE  PLACES  A  CURVATURE  ON  THE  WAVEFRONT  TO  BE 
C     USED  FOR  THE  CONVOLUTION  FORM  OF  FRESNEL  PROPAGATION. 
C 

C     CALL  QUAD1(NR,PI,DELMSH,DX,DY,WVL,DELX) 
C     ENDIF 
C 

1000  CONTINUE 
C 

C     THIS  SECTION  CALLS  OUT  GAUSSIAN  RANDOM  NUMBERS  FOR  THE  REAL  AND 
C     IMAGINARY  PHASE  ARRAYS. 
C 

CALL  GGAUS(NR,SEED) 
C 

C     THIS  BEGINS  THE  FILTERING  PROCESS  OF  THE  PHASE  SCREEN 
C 

CALL  FLTR(NR,DELMSH,CN2,DELX,WVL, FILTER) 
C 

C     HERE  THE  INDIRECT  TRANSFORM  IS  BEING  APPLIED  TO  THE  FILTERED 
C     PHASE  SCREENS. 
C 

CALL  IFTSCR(NR,M,DELMSH) 
C 

C     THE  FOLLOWING  ARRAYS  USED  IN  THE  FFT  ROUTINES  ARE  ZEROED  OUT  TO 
C     ENSURE  THAT  UNWANTED  VALUES  ARE  NOT  LEFT  IN  THE  ARRAYS. 
C 

DO  990  1=1, NR 
RE(I)=0.  0 
RIM(I)=0.  0 
990  CONTINUE 
C 

C     THIS  SECTION  DOES  THE  ALGEBRA  NEEDED  TO  MESH  THE  PHASE  SCREEN 
C     TOGETHER  WITH  THE  ELECTRIC  FIELD 
C 

DO  50  1=1, NR 
DO  50  J=1,NR 

XA=COS(PHASER( I , J) )*FIELDR( I , J) 
XB=COS(PHASER(I,J))-''FIELDI(I,J) 
XC=SIN( PHASER( I , J) )*FIELDR( I , J) 
XD=SIN(PHASER( I , J) ) *FIELDI( I , J) 
FIELDR(I,J)=XA-XD 
FIELDI(I,J)=XB+XC 
50  CONTINUE 
C 
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C     HERE  THE  FAST  FOURIER  TRANSFORM  IS  BEING  APPLIED  TO  THE  PERTURBED 
C     ELECTRIC  FIELD.  FOR  A  DIRECT  TRANSFORM  SIGN=-1.0,  AND  INDIRECT 
C     TRANSFORM  SIGN=+1.  0. 
C 

SIGN=-1.  0 

CALL  DFTIFT(NR,M,SIGN,DELMSH) 
C 

DO  141  1=1, NR 
DO  141  J=1,NR 
FIELDR(I,J)=FIELDR(I,J)/(NR*DELMSH) 
FIELDI( I , J)=FIELDI( I , J)/(NR*DELMSH) 
141  CONTINUE 
C 
C 

C     THIS  PORTION  OF  THE  IF  STATEMENT  CORRESPONDS  TO  THE  IMPLEMENTATION 
C     OF  THE  TRANSFER  FUNCTION  FORM  OF  FRESNEL  PROPAGATION.   THE 
C     SUBROUTINE  CALLED  TRNSFER  APPLIES  A  QUADRATIC  TO  THE  FIELD. 
C 

C     IF(NR.  LE.MODE)  THEN 

C     CALL  TRNSFR(NR,PI,DX,DY,WVL,DELX,DTRNS) 
C     ICNT=ICNT+1 
C     SIGN=+1.  0 

C     CALL  DFTIFT(NR,M,SIGN,DELMSH) 
C     GO  TO  888 
C     ELSE 
C 

C     THE  SUBROUTINE  CALLED  QUAD2  PUTS  THE  DIFFRACTION  PATTERN  IN  REAL 
C     SPACE  COORDINATES. 
C 

C     CALL  QUAD2(NR,PI,DX,DY,WVL,DELX) 
C     END IF 
C 

C  888  CONTINUE 
C 
C 

C     THIS  DO  LOOP  DETERMINES  THE  POWER  SPECTRAL  DENSITY  AND  SETS  IT  UP 
C     FOR  AN  FFT  TO  DETERMINE  THE  MCF. 
C 

DO  152  1=1, NR 
DO  152  J=1,NR 
FMCF(IJJ)=FIELDR(I,J)**2+FIELDI(I,J)**2 
FILL(I,J)=0.  0 
152  CONTINUE 
C 

C     ONCE  AGAIN  THE  ARRAYS  ARE  CLEARED  OF  STRAY  VALUES 
C 

DO  910  1=1, NR 

RE(I)=0.0 

RIM(I)=0.0 
910  CONTINUE 
C 

C     THE  INVERSE  FFT  IS  APPLIED  TO  THE  POWER  SPECTRAL  DENSITY 
C 

SIGN=+1.  0 

DO  901  1=1, NR 
DO  911  J=1,NR 
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RE(J)=FMCF(I,J) 
RIM(J)=FILL(I,J) 
911  CONTINUE 

CALL  FFT(M,SIGN,DELMSH) 
DO  921  J=1,NR 
FMCF(I,J)=RE(J) 
FILL(I,J)=RIM(J) 
921  CONTINUE 
901  CONTINUE 

DO  931  J=1,NR 
DO  941  1=1, NR 
RE(I)=FMCF(I,J) 
RIM(I)=FILL(I,J) 
941  CONTINUE 

CALL  FFT(M,SIGN,DELMSH) 
DO  951  1=1, NR 
FMCF(I,J)=RE(I) 
FILL(I,J)=RIM(I) 
951  CONTINUE 
931  CONTINUE 
C 
C 

C     THIS  SECTION  DETERMINES  THE  MAXIMUM  VALUE  AND  NORMALIZES  THE  MCF 
C 

FLDM=0.  0 
DO  89  1=1, NR 
DO  89  J=1,NR 
XMG=FMCF(I,J) 
IF(XMG. GT. FLDM)  THEN 
FLDM=XMG 
END  IF 
89  CONTINUE 
C     WRITE (*,*)  FLDM 
C     PAUSE 
C 

C     THE  MCF  IS  NORMALIZED  SO  THE  MAX  VALUE  IS  1.0 
C 

DO  29  1=1, N2 
DO  29  J=1,N2 
FMCF( I , J)=FMCF( I , J) /FLDM 
29  CONTINUE 
C 

C     THE  ATMOSPHERIC  MCF  IS  DETERMINED  BY  DIVIDING  OUT  THE  APERTURE 
C     FUNCTION  FROM  THE  PERTURBED  ELECTRIC  FIELD  MCF. 
C 

DO  91  1=1, N2 
DO  91  J=1,N2 
FMCF( I , J)=FMCF( I , J)/FNORM( I , J) 
91  CONTINUE 
C 

c 

C     DO  87  1=1,1 

C      DO  87  J=l,20 

C      WRITE(*,*)FMCF(I,J) 

C  87  CONTINUE 
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PAUSE 
C 

c 

C     THIS  SECTION  SETS  UP  THE  ARRAYS  TO  PLOT  THE  2-D  MCF 
C 

NT0T=2*NSIZE 
DO  827  1=1,1 
DO  827  J=l,NTOT 
DIST(J)=J 
FMAG(J)=FMCF(I,J) 
827  CONTINUE 
C 

C      DO  83  J=l,NTOT 
C       WRITE(*,*)FMAG(J),DIST(J) 
C   83  CONTINUE 
C     PAUSE 
C 

C     THE  FOLLOWING  SUBROUTINES  ARE  FOR  THE  GRAPHICS  PACKAGE 
C     THE  ATMOSPHERIC  MCF  IS  BEING  PLOTTED  AT  THIS  POINT 
C 

CALL  VSINIT(18,8.  ,10.  , 0 , ' MCF1. PLT' , IUNITV, IVID,5) 
CALL  ORIGINC  5,1.  5,0) 
CALL  SCALE(FMAG,5. ,NTOT,l) 

CALL  AXIS(0. ,0. ,'MCF' ,0,1,1,5. ,90. ,FMAG(NTOT+l) ,FMAG(NTOT+2) , . 1,1) 
CALL  SCALE(DIST,5. ,NTOT,l) 

CALL  AXIS(0.  ,0.  , '  ' ,0,-1,-1,5.  ,0.  ,DIST(NTOT+l) ,DIST(NTOT+2) , . 1,1) 
CALL  LINES(DIST,FMAG,NTOT,l,-l,ISYMB,. 1) 
CLOSE (IUNITV) 
C     CALL  INT86(ITNO,IREG) 

CALL  MSG(0.  ,0.  ,.15,' PRESS  ANY  KEY  TO  CONTINUE*, 0.  ,0,1) 
AN'S=GETC( ) 
CALL  GMODE(IVID) 
C 
C 

C     THE  MAGNITUDE  OF  THE  FFT'D  ELECTRIC  FIELD  IS  CALCULATED  IN  ORDER 
C     TO  PLOT  THE  OUTPUT. 
C 
C 

FMAX=0.  0 
DO  80  1=1, NR 
DO  80  J=1,NR 
FIELDM( I , J)=SQRT(FIELDR( I, J)**2+FIELDI(I, J)**2) 
X=FIELDM(I,J) 
IF(X.  GT.  FMAX)  THEN 

FMAX=X 
END  IF 
80  CONTINUE 
C 

C     THIS  SECTION  BEGINS  THE  CALLING  SEQUENCE  FOR  PLOTTING 
C 

CALL  VSINIT(MON,10.  ,8.  ,0, ' DITHER.  PLT' , IUNITV , IVID, 5) 
DO  100  1=1, N2 
DO  100  J=1,N2 
IX=J*INCR 
IY=I*INCR 
XMAX=ALOG10(FIELDM( I , J) /FMAX) 
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INDEX=8*ABS(XMAX)+1 

IF( INDEX. GE. 21)  THEN 
ICOLOR=0 
ELSE 

ICOLOR=NDEX( INDEX) 
END  IF 
CALL  PIXEL(IX,IY,ICOLOR) 
100  CONTINUE 

CALL  MSG(0.  ,1.  ,.  15, TITLE, 0.  ,0,0) 
CLOSE(IUNITV) 
C     CALL  INT86(1TN0,IREG) 

CALL  MSG(0. ,0. ,.15, 'PRESS  ANY  KEY  TO  CONTINUE', 0. ,0,0) 
ANS=GETC( ) 
CALL  GMODE(IVID) 
C 

C     THIS  IF  STATEMENT  QUES  THE  PROGRAM  TO  START  OVER  AGAIN  IF 
C     FRAUNHOFER  OR  THE  CONVOLUTION  FORM  OF  FRESNEL  ARE  USED 
C 

IF(NR. GT. MODE)  THEN 
GO  TO  999 
END  IF 
C 

C     THIS  IF  STATEMENT  QUES  THE  PROGRAM  TO  FULLY  COMPLETE  PROPAGATION 
C     THROUGH  ALL  THE  SLABS  IN  THE  TRANSFER  FUNCTION  FORM  OF  FRESENL 
C 

IF( ICNT.  NE. NUMSCR)  THEN 
GO  TO  1000 
END  IF 
GO  TO  999 
END 
C 

SUBROUTINE  QUAD 1 ( NR , PI , DELMSH , DX , DY , WVL , DELX) 

COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDI( 256 , 256) 

DX=DELMSH 

DY=DELMSH 

MID=(NR/2)+l 

DO  273  1=1, NR 

X=(I-(MID*2-1))*DX 
DO  273  J=1,NR 

Y=(J-(MID*2-1))*DY 

THETA=PI*( ( (X*X)+( Y*Y) )/(WVL*DELX) ) 
XX=FIELDR( I , J)*COS(THETA) 
YY=FIELDR(I,J)*SIN(THETA) 
ZZ=FIELDI(  I ,  J)-'-COS(THETA) 
VW=FIELDI(I,J)*SIN(THETA) 
FIELDR(I,J)=XX-WW 
FIELDI(I,J)=YY+ZZ 
273  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  QUAD2(NR,PI ,DX,DY, WVL, DELX) 

COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDI( 256 , 256) 

DX2=DX*WVL*DELX 

DY2=DY*WVL*DELX 

MID=NR/2+l 
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c 


DO  274  1=1, NR 
Y=(I-MID)--'DY2 
DO  274  J=1,NR 
X=(J-MID)'--DX2 

PHI=PI*(((X*X)+(Y*Y))/(WVL*DELX)) 
CX=FIELDR( I , J)*COS( PHI) 
CY=FIELDR(I,J)*SIN(PHI) 
CZ=FIELDI(I,J)*COS(PHI) 
CW=FIELDI(I,J)*SIN(PHI) 
CBR=CX-CW 
CBI=CY+CZ 

FIELDR( I , J)=(CBI/(WVL*DELX) ) 
FIELDI(I,J)=-1.*(CBR/(WVL*DELX)) 

274  CONTINUE 
RETURN 
END 

SUBROUTINE  TRNSFR( NR , PI , DX , DY , WVL , DELX , DTRNS ) 
COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDI( 256 , 256) 
DX=DTRNS/NR 
DY=DTRNS/NR 
MID=NR/2+l 
DO  275  1=1, NR 
FY=(I-MID)"DY 
DO  275  J=1,NR 
FX=(J-MID)"'DX 

FEE=-1.  *PI*WVL*DELX*(  (FX*FX)+(FY*FY) ) 

GX=FIELDR( I , J)*COS(FEE) 

GY=FIELDRCI,J)-SIN(FEE) 

GZ=FIELDI(I,J)»COS(FEE) 

GW=FIELDI(I,J)*SIN(FEE) 

FIELDR(I,J)=GX-GW 

FIELDI(I,J)=GY+GZ 

275  CONTINUE 
RETURN 
END 


SUBROUTINE  FLTR(NR,DELMSH,CN2 , DELX, WVL, FILTER) 

NOTE:  DELMSH'-'-NR  IS  THE  LARGEST  APERTURE  SIZE 

COMMON  /BLK3/  PHASER(256,256) ,PHASEI( 256 ,256) 

PI=3. 141592653589792 

POWER=-ll.  /6. 

TPI=2. *PI 

N2=NR/2 

NPIVOT=N2+l 

LAST=NPIVOT+l 

DLKAPA=( TPI/ ( NR*DELMSH) )**P0WER 

FACTOR=SQRT( (TPI**3)*.  033*CN2*DELX/(WVL**2) ) 

FUDGE=DLKAPA*FACTOR 

DO  100  I=l,NPIVOT 

EYE=REAL(I) 

EYE2=EYE*EYE 

DO    100   J=1,NR 
IF(J.  LE.NPIVOT)    THEN 
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WHY=REAL(J) 
ELSE 

WHY=REAL(NR-J+2) 
END  IF 
XKAPPA=SQRT( EYE2+WHY*WHY) 
AKAPPA=XKAPPA**FILTER 
PHASER ( I , J)=PHASER( I , J)*FUDGE*AKAPPA 
PHASEI( I , J)=PHASEI( I , J)*FUDGE*AKAPPA 
100   CONTINUE 

DO    110    I=LAST,    NR 
EYE=REAL(NR-I+2) 
EYE2=EYE*EYE 
DO   110   J=1,NR 
IF  (J.  LE.NPIVOT)   THEN 
WHY=REAL(J) 
ELSE 

WHY=REAL(NR-J+2) 
END  IF 
XK APP A=S QRT( E YE 2+WHY*WHY ) 
AKAPPA=XKAPPA**FILTER 
PHASER( I , J)=PHASER( I , J)*FUDGE*AKAPPA 
PHASEKI,  J)=PHASEI(I,  J)*FUDGE*AKAPPA 
110   CONTINUE 
C 

PHASER(1,1)=0.  0 
PHASEI(1,1)=0.  0 
C 

C     NOTE:  .  03441=TPI**-11.  /6.  TO  TRANSFORM  FROM  KAPPA  TO  FREQ  SPACE. 
DO  11  1=1, NR 
DO  11  J=1,NR 
PHASER(I,J)=PHASER(I,J)*.  03441 
PHASEKI, J)=PHASEI(I, J)*.  03441 
11  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  IFTSCR(NR,M,DELMSH) 
COMMON  /BLK1/  RE( 256) ,RIM( 256) 
COMMON  /BLK3/  PHASER(256 , 256) ,PHASEI( 256 ,256) 
PI=3.  141592653589792 
TPI=2. -PI 
SIGN=-1.  0 
DO  20  1=1, NR 
DO  21  J=1,NR 

RE(J)=PHASER(I,J) 
RIM(J)=PHASEI(I,J) 

21  CONTINUE 
CALL  FFT(M,SIGN,DELMSH) 
DO  22  J=1,NR 

PHASER(I,J)=RE(J) 
PHASEI(I,J)=RIM(J) 

22  CONTINUE 
20  CONTINUE 

DO  30  J=1,NR 
DO  31  1=1, NR 

RE(I)=PHASER(I,J) 
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RIM(I)=PHASEI(I,J) 

31  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 
DO  32  1=1, NR 

PHASER(I,J)=RE(I) 

PHASEI(I,J)=RIM(I) 

32  CONTINUE 
30  CONTINUE 

RETURN 
END 

SUBROUTINE  DFTIFT( NR , M , SIGN , DELMSH) 
COMMON  /BLK1/  RE( 256) ,RIM( 256) 
COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDI( 256 , 256) 
DO  60  1=1, NR 
DO  61  J=1,NR 

RE(J)=FIELDR(I,J) 

RIM(J)=FIELDI(I,J) 

61  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 
DO  62  J=1,NR 

FIELDR(I,J)=RE(J) 

FIELDI(I,J)=RIM(J) 

62  CONTINUE 
60  CONTINUE 

DO  70  J=1,NR 
DO  71  1=1, NR 

RE(I)=FIELDR(I,J) 
RIM(I)=FIELDI(I,J) 

71  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 
DO  72  1=1, NR 

FIELDR(I,J)=RE(I) 

FIELDI(I,J)=RIM(I) 

72  CONTINUE 
70  CONTINUE 

RETURN 
END 


SUBROUTINE  FFT(M, SIGN, DELMSH) 
COMMON  /BLK1/  RE( 256) ,RIM(256) 
PI=3.  141592653589792*SIGN 
N=2**M 
N1=N-1 
J=l 

DO  200  1=1, Nl 
IF(I.LT.  J)  THEN 

T=RE(J) 

RE(J)=RE(I) 

RE(I)=T 

T=RIM(J) 

RIM(J)=RIM(I) 

RIM(I)=T 
END  IF 
K=N/2 
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DO  201  WHILE(K.  LT.  J) 
J=J-K 
K=K/2 

201  CONTINUE 
J=J+K 

200  CONTINUE 
LE=1 

DO  202  L=1,M 
LE1=LE 
LE=LE+LE 
URE=1. 
UIM=0. 
ANG=PI/LE1 
WRE=COS(ANG) 
VIM=SIN(ANG) 
DO  203  J=1,LE1 
DO  204  I=J,N,LE 
IP=I+LE1 

TRE=RE(IP)*URE-RIM(IP)*UIM 
T I M=RE ( I P ) *U I M+R I M ( I P ) *URE 
RE(IP)=RE(I)-TRE 
RIM(IP)=RIM(I)-TIM 
RE(I)=RE(I)+TRE 
RIM(I)=RIM(I)+TIM 

204  CONTINUE 
T=URE*WRE-UIM*WIM 

U I M=URE * W I M+U I M*WRE 
URE=T 
203  CONTINUE 

202  CONTINUE 

IF(SIGN.  GT.  0.  0)  THEN 
PTS=1.  0/(N*DELMSH) 
DO  205  1=1, N 
RE(I)=RE(I)*PTS 
RIM(I)=RIM(I)*PTS 

205  CONTINUE 
END  IF 
RETURN 
END 


SUBROUTINE  GGAUS(NR , SEED) 

COMMON  /BLK3/  PHASER( 256 , 256) ,PHASEI( 256 , 256) 

DO  300  1=1, NR 
DO  300  J=1,NR 
301  V1=2.*RAN(SEED)-1 

V2=2.*RAN(SEED)-1 

S=V1*V1+V2*V2 

IF(S.GE.  1.  0)  GO  TO  301 

SCALE=SQRT(-2.*ALOG(S)/S) 

X1=V1*SCALE 

X2=V2*SCALE 

PHASER(I,J)=X1 

PHASEI(I,J;=0.  0 
300  CONTINUE 

RETURN 
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END 
C 

SUBROUTINE  MCF(FNORM,NR,M,PI ,NSIZE ,N2 ,DELMSH) 
COMMON  /BLK1/  RE( 256) ,RIM( 256) 
COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDI( 256 ,256) 
DIMENSION  FNORM(256,256) 
C 

C     THIS  SECTION  CREATES  THE  PLANAR  ELECTRIC  FIELD 
C 

DO  39  1=1, NR 
DO  39  J=1,NR 

FIELDR(I,J)=0.0 
FIELDI(I,J)=0.0 
39  CONTINUE 
C 

DO  45  I=N2-NSIZE+1,N2+NSIZE 
DO  45  J=N2-NSIZE+1,N2+NSIZE 
FIELDR(I,J)  =  1.  0 
45  CONTINUE 


SIGN=-1.  0 

CALL  DFTIFT( NR , M , S IGN , DELMSH ) 


DO  SO  1=1, NR 
DO  80  J=1,NR 
FNORM( I , J)=FIELDR( I, J)**2+FIELDI(I , J)**2 

80  CONTINUE 

DO  123  1=1, NR 
DO  123  J=1,NR 
FIELDI(I,J)=0.  0 
123  CONTINUE 

SIGN=+1.  0 
DO  90  1=1, NR 
DO  91  J=1,NR 

RE(J)=FNORM(I,J) 

RIM(J)=FIELDI(I,J) 

91  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 
DO  92  J=1,NR 

FNORM(I,J)=RE(J) 

FIELDI(I,J)=RIM(J) 

92  CONTINUE 
90  CONTINUE 

DO  93  J=1,NR 
DO  94  1=1, NR 
RE(I)=FNORM(I,J) 
RIM(I)=FIELDI(I,J) 
94  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 
DO  95  1=1, NR 
FNORM(I,J)=RE(I) 
FIELDI(I,J)=RIM(I) 
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95  CONTINUE 
93  CONTINUE 
C 

c 

FLDM=0.  0 

DO  88  1=1, NR 
DO  88  J=1,NR 
XMG=FNORM(I,J) 
IF(XMG.  GT.  FLDM)  THEN 
FLDM=XMG 

ENDIF 

88  CONTINUE 

DO  89  1=1, N2 
DO  89  J=1,N2 
FNORM(I,J)=FNORM(I,J)/FLDM 

89  CONTINUE 
C 

RETURN 
END 
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